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We report an essential improvement of the plain Fourier Monte Carlo algorithm that promises to 
be a powerful tool for investigating critical behavior in a large class of lattice models, in particular 
those containing microscopic or effective long-ranged interactions. On tuning the Monte Carlo 
acceptance rates separately for each wavevector, we are able to drastically reduce critical slowing 
down. We illustrate the resulting efficiency and unprecedented accuracy of our algorithm with a 
calculation of the universal elastic properties of crystalline membranes in the flat phase and derive 
a numerical estimate 77 = 0.795(10) for the critical exponent 77 that challenges those derived from 
other recent simulations. The large system sizes accessible to our present algorithm also allow to 
demonstrate that insufBciently taking into account corrections to scaling may severely hamper a 
finite size scaling analysis. This observation may also help to clarify the apparent disagreement of 
published numerical estimates of 77 in the existing literature. 
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Long-range interactions are ubiquitous in physics. Ex- 
amples include Coulomb, dipolar and higher multipole 
interactions in traditional condensed matter physics Jj, 
Wigner crystallization in fcrmionic quantum systems [2] , 
electrostatic interactions between cold trapped ions [3] 
or long-range elastic interactions between defects in con- 
densed matter systems [4 , just to name a few. Yet, even 
today many aspects related to long interaction ranges are 
only poorly understood [Sj- Sometimes, as e.g in many 
ionic systems, the long-range character is camouflaged by 
screening, leaving effective short-ranged interactions (see 
e.g. [5]). However, in a number of cases one is forced to 
deal with the full interaction range, frequently making 
theoretical attempts intractable and simulations compu- 
tationally expensive. 

Problems tighten further in simulations of critical long- 
range systems [71 [S] due to the required large system 
sizes and the notorious phenomenon of critical slowing 
down [9J. In recent years, cluster algorithms [TU] have 
been designed to overcome the latter problem. How- 
ever, they may be prohibitively difflcult to implement 
for complicated effective interactions arising e.g. in com- 
pressible spin models fllU12|. Yet, by utilizing the under- 
lying translation invariance, the structure of these effec- 
tive Hamiltonians often simplifies drastically on employ- 
ing the Fourier transform. Based on this observation, 
a radical approach was developed for lattice models in 
Refs. [HHS] and termed Fourier Monte Carlo (FMC) al- 
gorithm. Here we report a considerable improvement of 
this plain FMC algorithm, which practically also elimi- 
nates critical slowing down from the list of obstacles. The 
resulting optimized Fourier Monte Carlo (OFMA) algo- 
rithm is applied to study the elastic properties of solid 
membranes in the flat phase with unprecedented preci- 
sion. 

The present article is organized as follows. We begin 
with a discussion the manifestations of critical slowing 
down in plain FMC and an explanation of the optimized 



simulation scheme suggested by this analysis. This is fol- 
lowed by a short summary of the ideas underlying the de- 
scription of the asymptotic elastic behavior of fluctuating 
solid membranes and the related observables accessible 
in our simulations, and continued by a short explana- 
tion of our simulation setup. We present two approaches 
towards extracting a numerical estimate of the main crit- 
ical exponent r/ governing scaling behavior of solid mem- 
branes at long wavelengths from the generated data: (i) 
Investigation of the correlation function of out-of-plane 
deformations of the membrane gives only a preliminary 
estimate for 77 and reveals a peculiar finite size effect, (ii) 
A careful finite size scaling analysis of the membrane's 
mean squared displacement yields a presumably more re- 
liable numerical result and demonstrates the importance 
of properly taking into account subleading finite size cor- 
rections. The paper closes with a summary and short 
discussion of our results. 

OPTIMIZED FOURIER MONTE CARLO 

The original idea of Fourier Monte Carlo is quite sim- 
ple. In principle, any "spin" configuration {/(a;)} on a 
direct d-dim. lattice F, assumed to be real for simplicity, 
is in one-to-one correspondence with its set of (complex) 
Fourier amplitudes {/(g)} defined on the (first) Brillouin 
zone F. In FMC we completely forget about the direct 
lattice spins, treating the Fourier amplitudes f{q) as our 
basic Monte Carlo (MC) variables. A MC move consists 
of picking a random wave vector T and shifting 

/(g) ^ /(9) + e'^q.qo + ^*Sq,-q, , \t\ < T e (1) 

where e is randomly picked from a circle of fixed ra- 
dius centered around zero in the complex plane. The 
tricky part is, of course, how to compute the result- 
ing energy change IS.E accompanying the move ([T]) in 
an efficient way. Here we content ourselves with the 
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FIG. 1. Main plot; Dramatic reduction of integrated autocor- 
relation times rint(q) for moduli |/(<7)P for our new OFMC al- 
gorithm in comparison to plain FMC for system size L = 640 
and cutoff A = n/8. Deviations from linearity are due to the 
finite MC run time of 2^° MC steps (plain FMC) and finite 
tolerance ±5% for a 50% target acceptance (OFMC). Inset; 
g-dependent acceptance rates as measured for plain FMC. 



following brief description. Harmonic terms in a lat- 
tice Hamiltionian are diagonal under Fourier transform, 
and therefore it is straightforward to calculate the har- 
monic contribution to AE. In contrast, an anhar- 
monic contribution of type '^^/'^{x) turns into a sum 

Eqi...q4 Hqi) ■ ■ ■ /(94)Ar(qi + . . . 94), where the lattice 
delta function Ar{q) is defined to be 1 if q is a recipro- 
cal lattice vector and zero else. To avoid the resulting 
formidable combinatorial complexity, one trivially reor- 
ganizes J2xf^i^) = Sx(/^(^))^' which becomes diag- 
onal in terms of the Fourier amplitudes f^{q) of the 
squared field P{x). For a detailed account on how 
to efficiently calculate the anharmonic contribution to 
the energy change AE in terms of the amplitudes /(q) 
and P{q) and the nuts and bolts of FMC we refer to 
Refs. [TMS]. 

With FMC, the number of relevant degrees of freedom 
can be drastically reduced in calculating universal prop- 
erties if an effective Hamiltonian defined at wave vec- 
tor cutoff A is available. Consider e.g. a simple cubic 
lattice in d dimensions with lattice constant a. Then 
1 91 1 < Aq = 7r/a for any q G F. In this case the effort 
of an FMC simulation at cutoff A < Aq equals that of 
a direct lattice one with (Aq/A)"^ times more unit cells. 
Furthermore, as mentioned above, translation-invariant 
pair interactions are diagonal under a Fourier transform 
and thus pose no problem at all, regardless of their range 
- even Ewald summation is available (see Ref. |17)). 

Let us now turn to to explain the announced optimiza- 
tion of our algorithm. In the original version of FMC the 
radius for shifting the Fourier amplitudes in ([T]) is it- 
eratively optimized for an average acceptance rate (AR) 
of, say, 30 — 50% during the start-up of the simulation. 
Yet, even though the moves ([T]) are collective in nature, 
we observe a dramatic growth of the integrated auto- 
correlation times Tint(q) [Ml HH] of the squared ampli- 



tudes |/(q)P for modes close to the critical wave vector 
(taken to be qc = for simplicity). The sharp rise of 
Tintiq) for q=\q\^ 0, which is the crucial quantity that 
determines the statistical efficiency of measuring |/(q)P 
[T51 [TO] is the hallmark of critical slowing down. How 
can that be? The average amplitudes |/(q)| for q — > 
are much larger than those for ^ 0, but the algorithm 
attempts to move them all at the same maximum pace 
Tj. Thus the modes close to criticality simply make no 
headway in comparison to the noncritical ones, and one 
finds individual ARs close to 100% for the relatively few 
"small" q- vectors, while for larger wave vectors, the num- 
bers of which roughly increase as ^ q'^~^ , ARs drop to 
quite low values. However, nothing can prevent us from 
optimizing — r^{q) individually for each q in such a 
way that all modes /(q) separately enjoy the same uni- 
form AR. In practice, since changing r^[q) for one single 
q will influence all other individual ARs in a nonlinear 
way, we resort to a simple iterative procedure, aiming for 
a fixed collective rule-of-thumb target AR of, say, 50% 
with a tolerance of ±5% during the warm-up stage of the 
simulation. Excitingly, as soon as this initialization step 
is implemented, one observes an approximately uniform 
common value of T{at{q) with only weak q-dependence. 
Together with the collective nature of the move set ([I]), 
the dramatic suppression of critical slowing down makes 
this "Optimized FMC" (OFMC) algorithm an interesting 
alternative in cases where cluster algorithms are difficult 
to apply. 



SOLID MEMBRANES 

We illustrate the benefits of the abstract strategy out- 
lined above by considering the numerical determination 
of the exponent 77 governing the universal elastic prop- 
erties of solid membranes, a topic of high interest in its 
own right in molecular biology, medicine and pharmacy, 
chemical synthesis, and soft matter physics, just to name 
a few scientific disciplines. Owing to the recent meteoric 
rise of graphene this list has become even longer, 
including solid state physics, nanotechnology and elec- 
tronics. Due to space limitations, we make no attempt 
to do justice to all the sophisticated theoretical and com- 
putational approaches that have been developed to asses 
membrane elasticity and merely refer to the authorita- 
tive references [MI - EH] . Instead, we concentrate on the 
fact that for a crystalline membrane, which by defini- 
tion supports a nonzero static shear modulus /i 7^ 0, an 
effective long-range interaction between its out-of-plane 
deformations (OPDs) emerges as follows. 

Consider the so-called class of "phantom" membrane 
models for which self-avoidance effects are ignored. Liq- 
uid phantom membranes for which /i = 0, are known 
to collapse to a rotationally invariant "crumpled" phase 
characterized by an exponential decay of the membrane 
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unit normal correlations [H], and since only short-range 
interactions are at work, a transition to a "flat" phase 
via spontaneous breaking of this continuous symmetry 
is ruled out by the Mermin-Wagner-Hohenberg theo- 
rem |25j . For a crystalline membrane, however, elimi- 
nation of in-plane deformations (IPDs) from the parti- 
tion function by functional integration results in a shear- 
mediated effective long-range interaction between the 
OPDs which Mermin-Wagner-Hohenberg has nothing to 
say about [55]. And indeed, at sufficiently low tem- 
peratures crystalline membranes are found to be in a 
"flat" phase, in which the spatial correlations of the unit 
normals of the membrane tend towards a nonzero con- 
stant at long distances. In Refs. [27l |28], an effective 
Hamiltonian was formulated along these lines of think- 
ing. In the so-called Monge parametrization [21], defor- 
mations with respect to a given two-dimensional refer- 
ence plane with coordinates x = {xi,X2) are encoded 
in a height function f{x) parametrizing the OPDs and 
a two-dimensional vector u{x) of IPDs. Variations in 
J{x) give rise to a bending energy ^ J dPx (A/)^ [x] 
which is also present in liquid membranes, but for 
crystalline membranes f{x) also couples to the IPD's 
u(x) through an additional elastic stretching energy 
\ I d"^^ (2MAwf, (ic) -I- XAUii{x)ujj{x)) involving the 
Lagrangian strain tensor Uij{x) = [diUj{x) + djUi{x) + 
dif{x)djf{x)]/2. Transforming to reciprocal space and 
eliminating u from the partition function by Gaussian in- 
tegration [TT] [T^J HZl [5S] , one obtains an effective Hamil- 
tionian 

«.l/l^f/|f.,'l/ta)P + f/||l^(0)lt2) 

for the surviving OPD amplitudes f{q ^ 0), where 
Kfy — 4/iA(/iA + Aa)/ (2/iA + Aa) is the effective 2d Young 
modulus at cutoff A. The amplitudes J-{Q), which play 
a part similar to that of the amplitudes /^(q) of the 
squared field P{x) in the basic PMC algorithm outlined 
above, are defined via the nonlocal generalized convolu- 
tion in ^ 

HQ) ^= / X QfmnQ g) (3) 

whose wavevector dependence encodes the specific long- 
ranged interaction character (in ([s]) we have formally em- 
bedded the vectors Q and q in 3d). The basic observ- 
able for analyzing the thermodynamics resulting from 
(|2]), ([3| is certainly the correlation function of OPDs 
G{p)d^{p + q) = {f{p)f{q))- If anharmonic contribu- 
tions in ([2]) could be neglected, the equipartition theorem 
applied to the remaining harmonic bending contribution 
would yield G~^{q) = K\q^. In reality, the anharmonic- 
ity of ([2| causes the bare bending rigidity to be renor- 
malized, and ka picks up a nontrivial q-dependence, such 
that to leading order 

G-\q) = nt.{q)q\ ^A{q) ^ 9^" (4) 



The exponent 77 is the central quantity governing the uni- 
versal long distance elastic behavior of the flat phase. For 
instance, the mean square height fluctuations ((A/)^) di- 
verge for i — > 00 like 

((A/)2) = G(0) = j ^G(g) ^ L^C (5) 

with the roughness exponent [3D] C = 1 — 77/2, as can be 
seen from explicitly calculating the integral ([5| using an 
infrared cutoff \qi\ > 2tt/L and Q. 

SIMULATION SETUP 

In order to numerically determine rj, we performed 
OFMC simulations on a d = 2 square lattice and we mon- 
itored the set of squared moduli for wave vectors q 
inside a suitably chosen cutoff A, together with (A/)^ ^ 
J2q 1/(9) P the total energy E of the system. From 
the corresponding raw data series of these observables, 
which we also analyzed to ensure the complete equili- 
bration of our simulations before MC measurements, we 
calculated estimates of the corresponding statistical er- 
rors and the resulting integrated autocorrelation times 
T'mtiQ), T(Afy te using the jackknife approach [Hll9j. 
These estimates were cross-checked for consistency by di- 
rectly determining Ti„t('z), T{Af)2 and te from the raw 
data autocorrelation functions [TH (TS] and double-cross- 
checked using the blocking method [H |3T] . The choice of 
renormalized parameters ka = 0.1, K\ = 1.0, in which 
a factor {ksT)^^ has been absorbed, was motivated by 
the heuristic principle to have an approximate balance of 
harmonic and anharmonic contributions to total average 
energy changes. Compared to plain FMC, one indeed ob- 
serves the expected tremendous reduction of autocorrela- 
tion times Tint {q) (cf . Fig. fl| , T(a /)2 and te (not shown) . 
After completing this workwe realized that an optimiza- 
tion similar to ours had also been attempted in Ref . [32] , 
but was only implemented in the trivial case of a simple 
quasiharmonic model, and their collective "wave vector 
moves" were carried out at a frequency of one MC sweep 
on average, merely complementing a localized real-space 
MC move set. Suppression of critical slowing down with 
an efficiency comparable to that of our present method 
was therefore clearly out of reach. 

OUT-OF-PLANE DEFORMATION 
CORRELATIONS 

Unfortunately, trying to numerically extract rj from the 
correlation function G{q) of the OOP deformations via 
Q with sufficient precision suffers from several caveats. 
First of all, for G{q) no systematic FSS machinery seems 
to be available. Thus, even for a very large system size, 
simulation results may be contaminated by a residual 
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FIG. 2. Top: fit of simulation results for G{q) at linear system 
size L = 1800 and cutoff A = n/15. Error bars are smaller 
than symbol size. Bottom: Relative deviations of simulation 
data from the fit. Red lines are a guide to the eye. Gray lines 
indicate the relative statistical errors of the simulation data 
in percent. 



finite size dependence and effects of lattice anisotropy. 
Also, the quality of results may hinge crucially on the 
particular functional form of ansatz for the effective ex- 
ponent function required for interpolating between the 
mean field (MF) value rj = and the critical value of rj 
(cf e.g. Ref. ^). Integration of our ansatz 

77 



(6) 



1 + aq"^ 

with two parameters a and a besides rj to allow ad- 
justing both position and width of the crossover, yields 
the fit function G-^{q) = ^ ■ [l + {aq")-^]^'" . De- 
rived from a parallel version of our FMC code which 
will be described elsewhere. Fig. [2] shows the typical 
results of a corresponding fit obtained for a large sys- 
tem of linear size L = 1800 and cutoff A = 7r/15, ef- 
fectively resembling a direct lattice system of linear size 
Lo = 1800 X 15 = 27000. At first glance, the fit looks 
quite acceptable, producing a value of 77 = 0.751(17). 
However, comparing the relative statistical errors of the 
individual data to their relative deviations from the fit, 
our high statistical accuracy allows to resolve that most 
of the deviations remain within the range of the statis- 
tical noise of the data except for the smallest nonzero 
g-vectors (1, 0) • 2it/L and (0, 1) • 2'k/L on the cubic lat- 
tice, which exceed 10%. Following next in size are the 
vectors (1, 1) • 2ti/L and (1, —1) • 27r/L, who also show 
a noticeable deviation, albeit much smaller and oppo- 
site in tendency. To study the finite size dependence of 
these deviations, we decided to investigate a collection of 
systems with sizes L — 32rt, n — 1,2,..., 20, 22, . . . , 28 
and cubic cutoff A = 7r/8. Indeed, for growing system 
size L, the observed irregularities are qualitatively com- 
pletely similar but are shifted systematically in parallel 
towards limi_j.oo 2tt /L — (see Fig.|3]). Having ruled out 
trivial explanations for the observed behavior by various 
consistency checks, we conclude that what see is a fi- 
nite size effect related to the anisotropic structure of the 



convolution To extract a numerical estimate of rj 
from G{q), we thus decided to discard most deviatoric 
data, i.e. those for q \\ (1, 0), (0, 1), (1, 1), (1, -1). As 
Fig. |4] shows, the resulting fit is excellent and produces 
rj — 0.761(8). In view of the encountered difhculties, 
however, this estimate may be contaminated by resid- 
ual systematic errors beyond the pure statistical error of 
the fit we report here, as we shall argue below. Nev- 
ertheless, to appreciate the quality of our present data, 
note that the crossover ^ of G{q) from MF to criti- 
cal behavior roughly occurs at a "Ginzburg wave vector" 
P5] qc ~ y^3kBTKA/8nK\ analytically defined by the 
breakdown of the harmonic approximation. The ability 
of a simulation to efficiently sample the scaling region 
thus depends on the ratio p = qQ/{2TT/L) of qc to the 
smallest accessible wave vector component 2n/L. For 
instance, in the case of the atomistic MC simulations 
presented in Ref. |32| p turns out to be only approxi- 
mately 9. In contrast, for our largest systems we obtain 
p = qG/i2TT/L) w 66, thus providing p sa (66/9)2 _ 54 
times more data that actually explore the scaling region. 



OUT-OF-PLANE MEAN SQUARED 
DEFORMATIONS 



In contrast to G(q), analysis of ((A/)^) allows to em- 
ploy finite size scaling (FSS) techniques. On our discrete 
lattice, the integral in ([5| is replaced by a sum, whose 
asymptotic scaling behavior should comply to the gen- 
eral form By definition, an asymptotic scaling law 
of type ([5]) allows for various subleading algebraic and 
logarithmic corrections at finite L, we should therefore 
be included in a fit to the data in order to obtain pre- 
cise estimates for both rj and the corresponding error cr^ 
(see e.g. [S]). Unfortunately, however, to the author's 
best knowledge the structure of these corrections has not 
been worked out analytically up to date. We therefore 
have to allow for a priori unknown subleading corrections 
of logarithmic as well as power law type. In addition, an- 
alytical and numerical tests based on our crossover ansatz 
^ suggest the inclusion of a constant 6 > 0, such that 
we arrive at a FSS ansatz for (A/)^ ^ ^ 1/(9)1^ of type 



(A/)2 r^S + aL^-"{l + /3 In L + 7^-"') 



(7) 



in which we limit ourselves to including a single positive 
algebraic correction with exponent w. However, since 
the trade-off between logarithmic and algebraic correc- 
tions for small uj makes it numerically difficult to obtain 
meaningful fits, we content ourselves to studying both 
corrections separately. 

A purely logarithmic correction yields the value rj — 
0.793, but unfortunately this result does not inspire much 
confidence since the corresponding statistical fitting error 
is of the order of 10^. Turning to algebraic corrections, we 
fix 6 = and attempt to fit ([t]) to the data with a variable 
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correction exponent oj, which yields the vague result r] « 
0.781(100) for T], accompanied by the much too imprecise 
estimate a; = 0.372 ± 2.4 for A fit based on fixing w 
to this value produces r] = 0.784(5), but the dependence 
on the choice of uj remains unclear. To investigate this 
problem further, we thus decided to fit the data using 
the ansatz Q for a range of values of uj. We observe 
that only for roughly lo G [0.2, 1.2] these fits produced 
meaningful values for the parameters a, 7 and S within 
equally meaningful uncertainties. The results, which are 
gathered in Fig. [5] illustrate the unpleasant fact that 
the indeterminacy of u not only affects the quality of 
the error estimates for the sought-after exponent rj but 
indeed has a non-negligible effect on the estimated value 
of 1] itself. Indeed, Fig. [5] makes it obvious that without 
precise knowledge of the exponent uj, FSS based on the 
ansatz ([t]) does not allow to extract a reliable result for 
rj. Moreover, the lower part of Fig. (|6| demonstrates 
that, regardless of which value of w S [0.2, 1.2] is chosen, 
residual deviations between fits based on Eqn. ([T]) and 
the actual data at smaller L are found to persist. In 
combination with the pronounced w-dependence of the 
obtained values for ry, this provides compelling evidence 
that an ansatz of type ([t]) does not properly account for 
the finite size corrections to scaling. 

In contrast, an alternative fit based on the somewhat 
simple-minded ad-hoc ansatz 

{Af)^^6 + aL'-^{l + P/L + j/L^) (8) 

not only yields a comparably good agreement at large 
L (where any reasonable ansatz for (A/)^ with the cor- 
rect built-in asymptotics will works equally well) but ap- 
parently also captures the behavior at small L with an 
accuracy that seems very hard to improve any further 
(cf. again the lower part of Fig. ([6|). With this surpris- 
ingly good numerical agreement at both large and small 
L, it therefore appears to be of little relevance that the 
assumed representation of the scaling corrections as mere 
inverse integer powers of L lacks a strict theoretical jus- 
tification. Numerically, a fit using ([s]) produces the value 
r] — 0.795(5). Including a conservative safety margin of a 
factor of two in the error bar to this result, we are finally 
led to report the fair estimate 

r] = 0.795(10) (9) 



SUMMARY AND DISCUSSION 

In summary, we have explained how to optimize the 
plain Fourier MC algorithm and effectively eliminate crit- 
ical slowing down. As an application, we obtain high 
precision simulation results for the universal elastic be- 
havior of a crystalline membrane in the flat phase. For 
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FIG. 3. Comparison of relative deviations of simulation data 
for G{q) for various values of L and cubic cutoff A — n/8. 
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FIG. 4. Main plot: Same as top of Fig. [2] but with all contri- 
butions for q- vectors of symmetry type (1,0) and (1,1) omit- 
ted. Right upper inset: relative deviations of simulation data 
from the fit. Red lines are a guide to the eye. Gray lines 
indicate the relative statistical errors of the simulation data 
in percent. Left lower inset: crossover function mt. 



this problem, a considerable dispersion of previous es- 
timates of J] has been published in the existing litera- 
ture (cf. Table |l| for comparison with our present esti- 
mates) , and even a complete violation of scaling has been 
claimed [37]. While our present simulations fully sup- 
port the conventional universal scaling theory of solid 
membranes, our high precision and the accessibility of 
unprecedented effective system sizes allow to gain some 
new insight into the reason for the remaining numerical 
discrepancies in the published results for rj. Indeed, we 
find that the analysis of both G{k) and ((A/)^) requires 
an extremely careful analysis of finite size corrections. In 
the case of G(fe), finite size effects are serious enough to 
cast any attempt to extract rj from a naive scaling fit into 
doubt, but a finite size scaling analysis of ((A/)^) may 
also be severely biased by ignoring or incompletely taking 
into account subleading corrections. Our analysis reveals 
that such corrections clearly do affect the numerics of our 
present work even though the effective system sizes that 
we are able to access are much larger than those of previ- 
ous studies. We believe that these finding may also serve 
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TABLE I. Selection of numerical results for rj. Question marks 
indicate unreported error bars or possible systematic errors. 



V 



Method 



Rcf. 



0.2 0.4 0.6 O.i 



1.2 0.2 0.4 0.6 0.8 



0.750(5) MC, Gaussian spring pot., IPDs 

0.72(4) MC, Gaussian spring pot., OPDs 

0.849(?) nonperturbative RG 

0.821(?) self-consistent field approx. 

0.85(?) MC, atomistic carbon potential 

0.85(?) MC, MD, quasiharmonic model 

0.761(?) OFMC, G{q) 

0.795(10) OFMC, {(A/)^) 
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FIG. 5. Estimates and accompanying statistical errors ob- 
tained for the target exponent rj and the fit parameters a, 7, 5 
as obtained from fits of the FSS ansatz ^ for /3 = at various 
values of oj G [0.1, 1.2]. 
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FIG. 6. Top: fits of Eqs. ([7| for various fixed values om taken 
from the interval [0.1, 1.2]and Eqn. ([8| to simulation results 
for {{Aff) obtained for L = 32n, n = 1, 2, . . . , 20, 22, . . . , 28 
and A = 7r/8. Data error bars are smaller than symbol size. 
All fits are practically indistinguishable at this scale. Bottom: 
Relative deviations of simulation data from the fits. Lines are 
a guide to the eye. Gray lines indicate the relative statistical 
errors of the simulation data in percent. 



to explain the mentioned disagreement on the pubhshed 
numerical values for 77 in the existing literature that is 
immediately apparent from a glance at Table |lj 

The presented strategy for suppressing critical slowing 
is expected to work equally well for a large class of other 
lattice models at or near criticality. As to elastic mem- 
branes, the specific approach of the present paper should 
also be equally applicable to the hexatic case, for which 
we hope to obtain first results in the near future. 

We would like to thank K. Binder, C. Dellago, 
W. Janke, G. Kahl, M. Katsnelson, U. Pedersen and 
A. Travesset for discussions and acknowledge support by 
the Austrian Science Fund (FWF) Project P22087-N16. 
Major parts of our computations were performed on the 
Vienna Scientific Cluster (VSC2). 



* landreas.troe sterQtuw ien.ac.atI 
[1] N. W. Ashcroft and N. D. Mermin, Solid State Physics 

(Saunders College, Philadelphia, 1988). 
[2] A. Altland and B. D. Simons, Condensed Matter Field 

Theory (Cambridge University Press, Cambridge, UK, 

2010). 

[3] M. Lewenstein, A. Sanpera, and V. Ahufinger, Ultracold 
Atoms in Optical Lattices: Simulating quantum many- 
body systems (Oxford Unversity Press, Oxford, UK, 
2012). 

[4] P. M. Chaikin and T. Lubensky, Principles of Con- 
densed Matter Physics (Cambridge University Press, 
Cambridge, UK, 1995). 

[5] Long-Range Interacting Systems: Lecture Notes of the 
Les Houches Summer School, edited by T. Dauxois, S. 
Ruffo, and L. F. Cugliandolo (Oxford Unversity Press, 
New York, 2008), Vol. 90. 

[6] O. V. Patsahan and I. M. Mryglod, Journal of Physics: 
Condensed Matter 16, L235 (2004). 

[7] M. Fisher, Journal of Statistical Physics 75, 1 (1994), 
10.1007/BF02186278. 

[8] J. M. Caillol, D. Levesque, and J. J. Weis, Phys. Rev. 
Lett. 77, 4039 (1996). 

[9] D. Amit and V. Martin-Mayor, Field Theory, the Renor- 
malization Group, and Critical Phenomena: Graphs to 
Computers (3rd Edition) (World Scientific, Singapore, 
2005). 

[10] D. Landau and K. Binder, A Guide to Monte Carlo Sim- 
ulations in Statistical Physics, 3rd ed. (Cambridge Uni- 
versity Press, Cambridge, 2009). 

[11] D. Bergman and B. Halperin, Phys. Rev. B 13, 2145 
(1976). 

[12] A. Troster, Phys. Rev. Lett. 100, 140602 (2008). 

[13] A. Troster, Phys. Rev. B 76, 012402 (2007). 

[14] A. Troster and C. Dellago, Ferroelectrics 354, 225 (2007). 

[15] A. Troster, Comput. Phys. Comm. 179, 30 (2008). 

[16] A. Troster and C. Dellago, in The Fourier Monte Carlo 
Approach to Lattice Spin Models (Springer Verlag, Hei- 
delberg, Berlin, 2008), Chap, to appear in Computer 
Simulation Studies in Condensed Matter Physics XXI. 

[17] A. Troster, Phys. Rev. B 81, 012406 (2010). 

[18] K. Binder and D. W. Heermann, Monte Carlo Simula- 
tions in Statistical Physics, Graduate Texts in Physics, 5 
th ed. (Springer, Berlin, 2010). 

[19] B. A. Berg, Markov Chain Monte Carlo Simulations and 



7 



Their Statistical Analysis (World Scientific, Singapore, 
2004). 

[20] A. Gcim and K. Novosclov, Nature Mat. 6, 183 (2007). 
[21] Statistical Mechanics of Membranes and Surfaces, Second 

Edition, edited by D. Nelson, T. Piran, and S. Weinberg 

(World Scientific, Singapore, 1988). 
[22] M. Bowick and A. Travesset, Phys. Kept. 344, 255 

(2001). 

[23] M. 1. Katsnelson, Graphene: Carbon m Two Dimensions 
(Cambridge University Press, Cambridge, UK, 2012). 

[24] P. De Gennes and C. Taupin, J. Phys. Chem. 86, 2294 
(1982). 

[25] N. Mcrmin and H. Wagner, Phys. Rev. Lett. 17, 1133 

(1966). 

[26] L. Peliti and S. Lciblor, Phys. Rev. Lett. 54, 1690 (1985). 

[27] D. Nelson and L. Peliti, J. Phys. (Paris) 48, 1085 (1987). 

[28] D. Nelson, in Statistical Mechanics of Membranes and 
Surfaces, Second Edition, edited by D. Nelson, T. Pi- 
ran, and S. Weinberg (World Scientific, Singapore, 2004), 
Chap. 6. Theory of the Crumpling Transition, p. 131. 

[29] S. A. Safran, Statistical Thermodynamics of Surfaces, 
Interfaces, and Membranes (Perseus Books, Cambridge 



MA, USA, 2003). 

[30] M. J. Bowick, in Statistical Mechanics of Membranes and 
Surfaces, Second Edition, edited by D. Nelson, T. Pi- 
ran, and S. Weinberg (World Scientific, Singapore, 2004), 
Chap. 11. Fixed-connectivity membranes, p. 323. 

[31] H. Flyvbjerg and H. G. Petersen, The Journal of Chem- 
ical Physics 91, 461 (1989). 

[32] J. H. Los, M. I. Katsnelson, O. V. Yazyev, K. V. Za- 
kharchenko, and A. Fasohno, Phys. Rev. B 80, 121405(R) 
(2009). 

[33] E. Luijten and K. Binder, Phys. Rev. E 58, R4060 (1998). 
[34] M. J. Bowick, S. M. Catterall, M. Falcioni, G. Thorleif- 

sson, and K. N. Anagnostopoulos, J. Phys. I France 6, 

1321 (1996). 

[35] J.-P. Kownacki and D. Mouhanna, Phys. Rev. E 79, 
040101 (2009). 

[36] P. Le Doussal and L. Radzihovsky, Phys. Rev. Lett. 69, 

1209 (1992). 

[37] A. Fasolino, J. Los, and M. L Katsnelson, Nature Mate- 
rials 6, 858 (2007). 



